FAST OPED ALGORITHM FOR RECONSTRUCTION OF 
IMAGES FROM RADON DATA 



YUAN XU AND OLEG TISCHENKO 

Abstract. A fast implementation of the OPED algorithm, a reconstruction 
algorithm for Radon data introduced recently, is proposed and tested. The new 
implementation uses FFT for discrete sine transform and an interpolation step. 
The convergence of the fast implementation is proved under the condition that 
the function is mildly smooth. The numerical test shows that the accuracy of 
the OPED algorithm changes little when the fast implementation is used. 



1. Introduction 

Reconstruction images from Radon data is the central theme in medical imaging. 
The dominating reconstruction method, or algorithm, currently used for this task 
is FBP (filtered backprojection) method. A new reconstruction algorithm, called 
OPED, was introduced recently as a possible alternative in [12 and studied further 
in [9j[T3]. OPED is based on orthogonal polynomial expansions on the disk, which 
is fundamentally different from FBP. Our study has indicated a number of potential 
advantages of OPED vs FBP. The purpose of this paper is to show that OPED can 
compete with FBP in speed as well. 

Fast reconstruction is one of the reasons that FBP is widely used. The structure 
of the FBP algorithm contains a discrete convolution, which can be evaluated with 
FFT (fast Fourier transform). Together with an interpolation step, this ensures that 
FBP algorithm can be implemented efficiently. In contrast, the OPED algorithm 
does not use discrete convolutions. Nevertheless, it uses Chebyshev polynomials of 
the second kind, which turns out to allow an implementation with similar combina- 
tion of FFT and interpolation, except that we will use fast discrete sine transforms 
instead. The result is a fast implementation of OPED algorithm, comparable to 
the implementation of FBP in the number of operations. 

The additional interpolation step will introduce an error in the reconstruction. 
It turns out, however, that the additional error is small. The convergence of the 
OPED algorithm is proved theoretically under mild condition on the function that 
represents the image. Treating the approximation provided by the algorithm as 
an operator, denoted by Aim, from the space of continuous functions to itself, the 
proof in |12j amounts to show that the operator norm of ^ 2 m is exactly ||«42m|| = 
0(m\og(m+ 1)), where 2m+ 1 is the number of views of the Radon projections. It 
turns out, as we shall show, that the order of the norm remains to be the same even 
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when the interpolation step is used. The OPED algorithm preserves polynomials 
of degree 2m — 1, so that the order of ||^2m|| implies that the algorithm converges 
uniformly if the function is mildly smooth, say, if the function has continuous 
second order derivatives. The OPED with interpolation step no longer preserves 
polynomials. Nevertheless, we shall show that the convergence still holds for smooth 
functions, say for functions that have fourth continuous derivatives. We also test 
the fast implementation of OPED numerically. The result shows that the speed of 
fast implementation is improved in the order of magnitude, whiles the quality of 
the reconstruction changes little even for images that have sharp jumps. 

The paper is organized as follows. We describe OPED algorithm and its fast 
implementation in the following section. The convergence of the fast OPED algo- 
rithm is established in Section 3. The numerical testing and examples are given in 
Section 4. 



2. OPED ALGORITHM AND ITS FAST IMPLEMENTATION 

2.1. OPED algorithm. A two dimensional image on the unit disk B = {(x,y) : 
x 2 + y 2 < 1} is represented by a function f(x, y) defined on B. A Radon projection 
of / is a line integral, 

Kf{6,t):= I f(x,y)dxdy, < 9 < 2tt, -1 < t < 1, 

Jl(6,t) 

where 1(6, t) = {(x,y) : xcosO + ysm9 = t] n B is a line segment inside B. The 
essential problem of image reconstruction from Radon data is to recover the image 
represented by f(x, y) from a finite collection of its Radon projections. 

OPED is a reconstruction algorithm based on orthogonal expansion on the disk. 
Let V n (B) denote the space of orthogonal polynomials of degree n on B with 
respect to the Lebesgue measure. A function in L 2 (B) can be expanded in terms 
of orthogonal polynomials, that is, 

oo 

(2.1) /(z) = X> r °jfc/(z), Proj fe : L 2 (B) » V n (B). 

k=0 

The partial sum of this expansion is S n f(x) = X)fc=oP r °jfc f{ x iV)i which provides 
a natural approximation to /. It is shown in |12j that the partial sum can be 
expressed directly in terms of Radon projections: 

. 2m ^ -i 

(2.2) S 2m f(x 7 y) = ——Y / - nf{<f> v ,t)$ v {t;x,y)dt 
where <b - 2lxv 



2m+l ' 



^ 2m 

$ v (t;x,y) = — V7fc+ l)U k (t)U k (xcos(j) v + ysin^), 

Zm + 1 * — ' 

and Uk(t) denote the Chebyshev polynomial of the second kind, 
(2.3) U k (t) = Shlik+ a 1)d , t = cos9. 



It turns out that the formula (2.2 1 can be found implicitly in [3] (see (5.9), (4.3) 
and (3.7) there). Furthermore, a double integral expression for S n f can be found 
in [TJ|7] from which (2.2) can be deduced by a quadrature on the unit circle. 
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Using a quadrature to discretize the integral over t gives an approximation to / 
based on discrete Radon data. A particular choice of the quadrature formula gives 
the OPED algorithm: 

OPED Algorithm. Let m be a positive integer. For each reconstruction point 
(x,y) € B, 

2m 2m 

(2.4) Aqpbb ■= 'yi'Y] 9 j,vTj lV (x,y), g jtU =: Uf (0„,cos ipj) , 

i/=0 j=0 

where <j> v = ipj = and 

^ 2m 

(2.5) Tj, v (x,y) = — — — VVfc + l)sin(fc + jU k {x cos <f)„ + ysin^). 

(2m + l) 2 f^ 



The above version of OPED with be refereed to as OPED of type I as it comes 
from Gaussian quadrature based on the zeros cos ^jrjp ; < j < 2m, of the 
Chebyshev polynomial of the first kind. An alternative is the OPED of type II, 
which uses Gaussian quadrature based on the zeros, cos „ 3 * , 1 < j < 2m, of 



the Chebyshev polynomials of the second kind, for which the sum over j in (2.4) 



should start from j = 1 and ipj should be replaced by ipj = 2m \i • These two 
choices lead to different scanning geometry of x-ray data, the advantage of the first 
type is explained in [S]. Unless specifically stated, OPED in the following means 
type I. The type II case will be mentioned whenever appropriate, the proof can be 
worked out in the same way and will be omitted. 

For each (x, y) in the domain B, the algorithm produces a number .Aoped which 
approximates the value of the image at that point. Typically a larger m produces 
a better approximation. The values of -4oped on the reconstruction points, say on 
a grid of pixel points, give the reconstruction of the image. The convergence of the 
OPED algorithm is equivalent to the approximation property of the linear operator 

2m 2m 

(2.6) A 2m f(x,y) = Y,J2 n f^ 

cos ipj)Tj :U (x,y) 

where / is the function that represents the image. The operator preserves 
polynomials of degree 2m — 1, that is, A2mf = / if / is a polynomial of degree at 
most 2m — 1. In other words, the algorithm produces an image exactly if the image 
happens to be represented by a polynomial of degree less than 2m— 1. Furthermore, 
let ||^2m|| denote the operator norm of A.2m hi the uniform norm over B; then f|12p 

(2.7) \\-A 2m \\ =0(mlog(m + l)), as m -> oo. 

As a consequence of this estimate, it follows that A 2m f converges uniformly to / 
on the disk B if / has continuous second order derivatives. The numerical test 
has shown that OPED gives fairly accurate reconstruction even when the data has 



sharp singularities. The proof of (2.7) is based on the following compact formula 
of T jiV (x,y). 
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Proposition 2.1. Let 6 V := 9 u {x,y). Then 

(0 ( ^ - (- 1 ) J (2m + 1) sin(2m + l)fl„ (-l)^(2m + 1) cos(2m + 1) 

[ + j J> ^' yj ~ 2 sin ^ ~ cos ^ - cos 0, 

sin sin ^ — (— 1) J sin(2m + 1)^(1 — cos 0j, cos i/>j) 
2sm0 u (cosipj — cos6 u ) 2 



It should be mentioned that (2.7 1 was proved in [T^] for OPED of type II. The 
compact formula of Tj tV in [12] is also established for type II. So, the formula of 
Tj u in the above proposition and (2.7 1 are in fact new, but the proof is very much 



similar to that of [T^] and we choose not to repeat it. 

The use of orthogonal expansion in reconstructing functions from their Radon 
data can be traced back to the classical paper of Cormack, whose Radon inversion 
formula is based on the expansions of / and Rf in spherical harmonics ([6j p. 25]). 
The orthogonal polynomials on the disk have been used in the papers [3j [5] and in 



the singular value decomposition of Radon transforms (see [B]). The formula (2.2) 



and the resulted OPED algorithm are introduced in [T5]. As we mentioned before 



that ((Z2J) itself can be deduced from (3j as well as BG,P. The OPED of type II 
is closely related to an algorithm in f2j. The connection to orthogonal polynomial 
expansion was not considered in [2], nor was the convergence studied there. For 
other algorithms and image reconstruction based on polynomial approximation, see 
[21 SI 03 E] and the references therein. 

2.2. Fast implementation of OPED algorithm. The structure of the A^m 



in (2.4) allows us to use FFT (fast Fourier transform) once in a straightforward 



manner. In fact, let us define 

k + 1 2m 
Sk ' u = (2m + I) 2 ^ 9j ' u Sin(/C + ' 

Then Sk, v can be evaluated by FFT for discrete sine transform. We can write 
-^oped as 

2m 2m 

(2.8) -4oPED = ^$> fc ,vUk (% cos (/>„ + y sin fpy) . 

Thus, the main step of the OPED algorithm lies in the evaluation of the above 
double sum, which can be considered as a back projection step. Let N = 2m + 1. 
Then the evaluation of the matrix Sk,u costs 0(N 2 log N) operations (flops). The 
evaluation of the double sum costs 0(N 2 ) operations. Hence, the cost of evaluation 
on a grid of M x M is 0(N 2 (M 2 + log N)). In particular, if M w N, then the cost 
is 0(N 4 ). The main operation cost is at the evaluation of the double sums at the 
grid points. In other words, the main cost lies in the back projection step. 



Unlike the FBP algorithm, the sum in (2.8) does not contain a discrete convo- 
lution that can be evaluated via FFT. However, the formula of Uk in (2.3) allows 
us to write 

2m _^ 2m 

(2.9) A peb(x, y) = J2 ~YT ^ Sk ' v + 

v=0 v fc=0 

where 

6 V := 6 v (x,y) = arccos(i£ cos (f> v + ?/ sin <?!>„). 
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The inner sum can be evaluated by FFT for discrete sine transform at certain 
points, which suggests that we introduce an interpolation step to take advantage 
of the fast evaluation by FFT. To be more precise, we define 

2m 
k=Q 



after the inner sum of (2.9 1. The FFT for discrete sine transforms can be used to 
evaluate the numbers 

oli, v := a„(6), & := ^m+l 1 * = 0, 1, . . . , 2m - 1 



effectively. That is, the inner sum in (2.9 1 can be evaluated by FFT when 9 v (x, y) — 
fj. For the interpolation step, we choose linear interpolation, which is also used 
in the implementation of FBP [6, p. 109]. For a given (x,y), we choose the 
integer I such that Q v (x, y) lies between 6 and 6+i an d use the value of the linear 
interpolation between ai tV and a/+i w as an approximation to the inner sum of (2.9 1. 
The linear interpolation is given by 

(2.10) £ u (9)=u u (9)ai +1 , u + (l-u v (9))ai, 1/ , u„(fl) := / ~ , 

6+i - ?/ 

where 6 < 9 < 6+i- Then the fast implementation of OPED is given as follows: 
Fast OPED algorithm: Let m be a positive integer. 

Step 1. For each v = 0, . . . , 2m, use FFT to compute for each k = 0, . . . , 2m, 

k- 



^^>sin(fe + l)^ 



(2m + l) 2 ^ 



2. For each I = 0,1, ... , 2m — 1, use FFT to compute 

2m 

OLi. v := Sk, v sin(fc + 1)6- 

fc=0 

Step 3. For each reconstruction point (x,y) inside the disk of the radius cos 2 m+i > 
determine integers I such that 

9 U —1, where 9 V = arccos(a;cos § v + sin^), 



I = 

and evaluate 



2m ^ 

/oped = ^ . „ [(1 - u v )ai, v + u v ai + x , u ] , 
^-^ sm 9 V 

where u v = (2m + l)9 v /ir — (I + 1). 

A fast implementation for OPED of type II works similarly, with 6 replaced by 
£ = ^2»n+T) < ^ < 2m. In other words, exchanging the values of ipj and 6 in the 
above fast algorithm for OPED of type I leads to the fast algorithm for OPED of 
type II. 

A couple of remarks are in order. First of all, sin 9 V appears in the denominator in 
the last step of the algorithm. However, sin 9 V = only if cos 9 V (x, y) = x cos <j) v + 
ysin0„ = 1, which happens only if (x,y) — (cos <j> v , sin <j> v ) . Since the points 
(cos</>„, sin (f> v ) are on the boundary of the region B, we do not have to take them 
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as reconstruction points. In fact, the region of interests is usually inside the unit 
disk; thus, we can evaluate at points inside a smaller disk in B, and this will also 
ensure that the values of sin B v in the last step will not be too small to cause loss 
of significant digits in the computation. Furthermore, if we restrict (x, y) to a disk 
with radius cos£o = cos7r/(2m+ 1), then it will also guarantee that the choice of I 
in the Step 3 is unique for all (x, y) in that disk. 

The algorithm uses FFT twice, the final sum in Step 3 is a single sum whose 
evaluation costs O(N) operations. Hence, the cost of evaluations on an M x M grid 
with M k, N is 0(N 3 ), which is in the same order of magnitude as the FBP algo- 
rithm. In fact, the structure of the OPED algorithm with this fast implementation 
is similar to that of FBP algorithm (cf. [6l p. 109]). 



3. Convergence of OPED algorithm with linear Interpolation 

Throughout this section, we let || • | denote the uniform norm on B. 

As we mentioned before that the convergence of the OPED algorithm depends 
on the norm A 2m of the operator defined in (2.6). Following the proof in [T2] . we 
have 

2m 2m 

1 1 -42m 1 1 = max V" V" sin tpj \T jM (x,y) | = C(mlog(m+ 1)). 

As a consequence of this estimate and the fact that A2 m preserves polynomials of 
degree 2m — 1 , the triangle inequality implies that 

(3.1) \\A 2m f-f\\ <cmlog(m+l)E 2m -i{f) 

where E n (f) = inf{||/ — P n \ : degP n < n} is the error of best approximation to / 
by polynomials of degree at most n. It is known ([H]) that if / £ C 2r (B), then 

(3.2) E n (f) <cn- 2r \\V r f\\, 



where V is a differential operator of order 2r. In particular, (3.1 I and (3.2 1 show 
that A2 m f converges uniformly to / if / has continuous second order derivatives. 

When the interpolation step is used in the evaluation of OPED, the operator is 
changed and we denote the new operator by Al2m, which is given by 

2m ^ 

Al 2m f(x,y) := -■ (n ( y, 4 J/)), 

,,-n sll H t 'iH a '' V)) 



where l v is defined in (2.10). It turns out that the norm of the operator AT 2m has 
the same growth order as that of ^m- 

Proposition 3.1. For m > 0, let £! TO = {(x, y) : yj x 2 + y 2 < cos7r/(2m + 1)}. 
Then 

max \AT 2m f(x,y)\ < c||/|| m(log(m+ 1)). 

Proof. Let u v = (& v — — The way that the index I is chosen implies 

that < u v < 1. Recall the definition of Tj^(x,y) in (2.51. We shall abuse the 
notation somewhat and write Tj <v {^i) when arccos(xcos 4> w + ysm(f>v) — ^. Using 
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the definition of U n (t) in (2.3 1, the operator can be written as 

2m ^ 2m 

i/=0 v ^ ' w ^ j=o 

x [(1 - u v )Tj >v (£i) sin. & + u I/ T J - iI/ (^ +1 )sin^+i] 



For (a;,y) G f2 TO , write x — rcos<f> and y — rsxncf), then i cos ^„ + ysin( 
rcos(0 — <j) v ) < r and r < cos7r/(2m + 1). In particular, 



(3.3) 



sin 6„(x, y) = y/l - cos 2 6 v {x,y) > Jl - cos 2 = sin . 



Furthermore, the definition of shows that there is a constant c independent of m 
such that 



sin & 



< c 



and 



sin 6 



z+i 



sin(0„(:c,y)) 

for (x,y) g fi m . Consequently, we conclude that 



sia(6 v (x,y)) 



< c 



(3.4) \M 2m f(x,y)\ <J2J2\ n f(^ C0S ^)\ (\ T iAti)\ + l^>(& + i)|) ■ 

!/=0 j=0 

Using the fact that \TZf(<f>,t)\ < yl — i 2 ||/|| (see, for example, [H]), it follows that 



\AJ 2m f(x, y )\ < n/ll J2J2 sin ^ (\ T iAZi)\ + |r J >te+i)|) 

i/=0 3=0 



We now use the compact formula of 7} „ in Proposition 2.1 Since sin(2m+ 1)£; = 
and cos(2rn+ 1)£ ; = (-l) z+1 , it follows that 



(2m+l) 2 T^(6) 



(-l)^+i+ 1 (2m + l) (-l) /+1 sin^- 



COS Ipj — COS £z 2(C0S ijlj — COS £;) 2 



Since 
have 



_ (j+l/2)x 



} 2m-tr an d £z = 2m+T ' * ne denominator of !}.„(£;) is never zero. We 
cos i/jj — cos = 2 sin — - — - sin — - — - . 



Since ip2m-j = tt — ipj and ^m-i = n — we can assume that < < 7r/2, 
which means that < I < m — 1. If < V'j < 7r /2, we have 



sin sin % cos tf ib n 

3 =2 2 , 2 < 2 cos < 2. 



The above equation also holds if n/2 < ipj < n, since then 7r/4 < ' < 37r/4 
and sin > \/2/2. Hence, using the fact that sinfl > (2/tt)0 for < 9 < tt/2 
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we have 

2m 2m 2m 



Z^Z_^ yji j.mwi - Z^ | ■ V- 3 -6 i 2m + l^ 



^=o j=o ,-=o I sm 2" I 1 i=o I sm — 



2 



2m o 2m ^ 

" 1^-61 + 2m+l^ 

^\ 2m + 1 ^\ 2m + 1 
^|j_Z+l/2| + ^|j-Z+l/2|2 

< c(2m + l)log(m + 1). 



The sum involving Tj_ v {^i + \) is estimated in exactly the same way. By ( 3.4 ) , the 
proof is completed. □ 

This proposition shows that the interpolation step does not increase the growth 
order of the operator norm, at least when we restrict the norm to the region £l m . 
Unlike A^m, the operator Al2m no longer preserves polynomials of high degrees. 
We need an estimate of the error AX-2mf — f for / being a polynomial. First, 
however, we need some results on proj fe : L 2 (B) i— ► Vk{B) defined in (2.1). Let 
P„(x,y), x,y £ B, denote the reproducing kernel of Vfc(P). Then 

proj fc /(x) = - / /(y)P fc (x,y)dy. 

Furthermore, the kernel Pj. satisfies a compact formula ([TU]) 

(3.5) P fe (x,y) = *±1 £ U k ((x,y) + s^WV^WF) , 

where (•, •) and |j • || are the usual Euclidean inner product and norm on M 2 , and 
Uk is the Chebyshev polynomial of the second kind. 

Lemma 3.2. Let \\fW2 denote the L 2 {B) norm of f. Then the projection operator 
satisfies 

|proj fe /(x)|<(fc+l)||/|| 2 , VxeP. 

Proof. Using the Cauchy-Schwartz inequality and the reproducing property of Pfe(-, •), 
we have 

1/2 



|proj fe /(x)| < ||/|| 2 ^^|P fe (x, y )| 2 dy) = ||/|| 2 [P fc (x,x)] 



1/2 



The definition ( |2.3[ ) gives the well known inequality |L/fc(t)| < k+1, so that by (3.5) 
we have |Pfc(x, x)| < (k + l) 2 , from which the stated result follows. □ 

Using proj fc / we can construct a sequence of polynomials that approximates /. 
For this we let 77 be a nonnegative C°° function on E that satisfies 

*?(*) = 1, 0<i<l, and supp7?(*) = [0, 2]. 

For / e L 2 (B) we define S%f by 

2n / ■ \ 

SZf(x) = J2v[^)proi j f(x). 
3=0 ^ n ' / 
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Let LT^ denote the space of polynomials of degree at most n in two variables. The 
approximation property of S% is given in the following lemma ([11 ). 

Lemma 3.3. Let f e C(B). Then 

(1) S%f e nL_i and S%P = PforPe U 2 n ; 

(2) for n E N, 

ll^/H < c||/|| and \\f-SZf\\<cE n (f). 

This shows that S%f realizes, up to a constant, the best approximation by poly- 
nomials. For this sequence of polynomials, we have the following: 

Lemma 3.4. Let F = Sfj and q < m. For each compact subset £1 in B, there is a 
constant c, independent of q, m and f , such that 

max \AI 2m F(x,y) - F(x,y)\ < c / — — 

Proof. Since F is a polynomial of degree 2q, we write its expansion as 

F(x, y) = J2 F i( x > V)> F i = r )\- a ) P r %" /(*) e V ^ B )- 

3=0 

The Radon projections of orthogonal polynomials can be explicitly computed. We 
have ([5]), 

nFj(<t>,t) = -r^py s/l - t 2 Uj (t) Fj (cos <f>, sin <j>) . 

The polynomials Uk(t) are orthogonal with respect to yl — t 2 on the interval 
[—1,1]. Consequently, 

If 1 .„ . 1 



KF 3 (4>,t)U k (t)dt= — — i^(cos0,sin<£)(S fc j. 

7T J ~r t 



Hence, using (2.2) and the fact that q < m, we conclude that 

(3.6) F(x, y) =S 2m F(x, y) = S 2m Fj (x, y) 

3=0 

^ 2m 2<? 

= 2m+ i E E F J ( cos 0" ' sin <^ ) ^ (cos 0„ (a;, y)) . 
i/=o j=o 

In exactly the same way, we also obtain that 

2m ^ 2q 

(3.7) ^/(a, y) = ^—^ £ g ( } E *i ( C0S ^. Sin *0 

I/=0 i J=0 

x [(1 - u„) sin(j + l)^ + i + u v sin(j + 



where u u = u v {9 u {x 1 y)) is the linear function in (2.10). Introducing the function 

2q 

sin</v) sin(j + 1)0, 
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which depends on F, we then conclude that 
F(x,y) - Al2 m F(x,y) 

^ 2m ^ 

= o — XT E —m \ [Hv(0v(x, V)) - (1 - Uu)H u (fy) - u v H v {i l+1 )] . 

£Tn -\- i ^ ^ sin t/jy^x, y j 

The expression in the square bracket is the difference of H v {-) and its linear inter- 
polation at and evaluated at 9 v {x, y). Hence the well-known estimate of the 
error of linear interpolation shows 

[H„{6 v {x, y)) - (1 - u v )H v {^) - u v H v ^ l+x )} 

< h\HZ\\\(e v (x, y )-mWx,v)-M\ < ilKII 



2 ,i v \uv»\~,»j _ 16 „ ^ii ( 2m + 1)2' 

where = maxo^e^Tr |iJ"(6>)|. Since ifi, is a trigonometric polynomial of degree 

2q, the classical Bernstein inequality (cf. [TH Vol. 2, p. 11]) shows that 

\\H'J\\ < (2q) 2 \\H u \\ :=2q 2 max \H u (t)\. 



0<t<1T 



Furthermore, by the definition of H v and Lemma 3.2 we have 

2q 2q 

llH^K^lFjix^lKj^lproijF&y)] 

3=0 j=0 

<ll/l| 2 E(j + l)<2(g+l) 2 ||/| 
j=o 

Hence, putting these inequalities together, we conclude that 

4 2m 1 

\F(x,y) -AL 2m F(x,y)\ < c\\f\\ J ^ 

(2m + 1 r 



(2m + 1) 3 ^sin^^.y) 



If (a;, y) £ f2, a compact set in B, then there is an r < 1 such that y/ x 2 + y 2 < r, so 
that svn.9 v (x,y) > %/l — r 2 as in (3.31. Thus the last sum is bounded by a constant 
multiple of m and the stated result follows. □ 

Theorem 3.5. If f £ C 4 (B) and Q is a compact subset of B, then there is a 
constant Cf, independent of m, such that 

(3.8) max \AI 2m f(x, y) - f(x, y) \ < c f l ° g{m + 1} . 

(x,y)en V™ 

In particular, AI 2m f converges uniformly to f on any compact subset of B. 
Proof. Let F = S^j f and q < m as in the previous lemma and be a compact subset 



of B. Since AX2m is a linear operator, by the triangle inequality, Proposition 3.1 
Lemma 13.31 and Lemma 13.41 we obtain 



\AX 2m f(x,y) - f(x,y)\ < (1 + ||^4X 2T „||)j(/- J F)(ar,y)| + \Al 2m F(x,y) - F(x,y)\ 

< c (mlog(m + l)E q {f) + ^- 2 g 4 ||/||) 
for all (x,y) S f2 and q < m. In particular, if / G C 4 (B), then (3.2) implies that 
(3.9) max \Jtt2 m f(x,y)-f(x,y)\<c(mlog(m + l)q- 4 \\V 2 f\\+ m - 2 q 4 \\f\\). 



In particular, setting q w m 3 / 8 in the above inequality gives (3.8 1. □ 
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The theorem is stated for functions in C 4 (B), which is likely too restrictive and 
the convergence could hold under less restrictive conditions. Also, for functions 



that are more smooth, one could get better convergence rate by choosing q in (3.9) 



differently. The rate so obtained and the one in (3.8 1, however, are likely not sharp 



and reflect the limitation of our method of proof. 

The main merit of the theorem is that it establishes the convergence of OPED 
algorithm with linear interpolation for smooth functions (images). In practice, 
the images often have sharp edges, which means that the functions representing 
images may not be smooth or even continuous. In the following section we present 
numerical examples, which demonstrate that Al^mf, OPED with interpolation 
step, converges well even for functions that are not continuous and it converges 
almost as well as A2 m f, OPED without interpolation step. 

4. Implementation and Result 

For numerical implementation, we used the FFT for discrete sine transform in 
the package FFTW ( |http:/ /www . fftw.org/ ). The numerical example is conducted 
on the Shepp-Logan head phantom [5] (see Figure 2). This is an analytic phantom, 
highly singular, as the image contains jumps at the boundary of every ellipse in the 
image, include the one on the boundary. The function that represents the image is 
a step function, which is not continuous at the boundary of each of the ellipses. 

We reconstruct the image with OPED algorithm without the interpolation step 
and Fast OPED algorithm, which contains the interpolation step as shown in pre- 
vious section, respectively. In both cases, Sk,u are computed with FFT. 

In Figure 1 images reconstructed by the original OPED and Fast OPED algo- 
rithms are depicted side by side. In these images we choose m — 512 and the size 
of the images are 512 x 512 pixels. 




Figure 1. Left: reconstruction by OPED algorithm with m = 512. Right: recon- 
struction by Fast OPED algorithm with m = 512. 

The reconstruction is carried out on a CELSIUS R610 computer with two In- 
tel Xeon(TM) CPU, each 3065 MHz, and 4 GB RAM. The code is written in C 
language. Using OPED algorithm, the reconstruction took 344 seconds, in which 
more than 95% of the time is used on the back projection step. Using Fast OPED 
algorithm, the reconstruction took merely 13 seconds, an improvement of more that 
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26 times. Furthermore, the two images show almost no visual difference. In Figure 
2, the original Shepp-Logan phantom and the difference of the two images in Figure 
1 are depicted. 



Figure 2. Left: original Shepp-Logan phantom. Right: the difference between the 
two reconstructed images in Figure 1. 

We can also measure the errors of the reconstruction. Let X denote an image, 
represented by its pixel values, so that X = {JQ : 1 < i < N}, where N is the 
number of pixels in the image. Let X R denote the reconstructed image, X R = 
{X R : 1 < i < N}. The relative square total error (RSE) between X and X R is 
defined by 



In our case N — 512 x 512 = 512 2 . The result is reported in Table 1, in which ORIG 
stands for the original Shepp-Logan phantom, OPED and Fast OPED stand for 
reconstructed image by OPED and by Fast OPED, respectively. For example, ORIG 
vs OPED means the error between the original phantom and the reconstruction by 



and the mean error (ME) between X and X R is defined by 





OPED. 



ORIG vs OPED orig vs Fast OPED OPED vs Fast OPED 



RSE 



0.00239702 0.00249574 0.000515499 



ME 



0.0129175 0.00981329 0.007715128 

Table 1. Error Estimates of OPED and Fast OPED 



The results in the table shows that fast OPED is slightly worse in relative least 
square error, but slightly better in mean error. The order of magnitude of the error 
is the same. The difference is practically negligible. 
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5. Conclusion 

We introduced a fast implementation of OPED algorithm by using FFT and a 
linear interpolation step. The fast algorithm is proved to converge uniformly on any 
compact subset in the disk and the numerical test has shown that it reconstructs 
images accurately and is as good as the original OPED algorithm. 

In conclusion, the fast OPED algorithm with interpolation step is not only much 
faster, it also shares main merits of the original OPED algorithm. Hence, for prac- 
tical applications, the fast OPED algorithm should be the recommended method. 
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